Dynamic  Models  Including  Uncertainty 


AFOSR  FINAL  TECHNICAL  REPORT 
FA9550-08-1-0147 

For  the  period:  February  1,  2008  November  30,  2008 

submitted  to 


Air  Force  Office  of  Scientific  Research 
Attn:  Dr.  F.  Fahroo  and  Dr.  W.  McEneaney,  Dynamics  and  Control 
Phone  (703)  696-7796  Fax  (703)  696-8450 

875  N.  Randolph  St,  Suite  325,  Room  4111 
Arlington,  VA  22203-1768 

January  22,  2009 


PI: 

Dr.  H.T.  Ranks 

North  Carolina  State  University 
Center  for  Research  in  Scientific  Computation 
Box  8212 

Raleigh,  NC  27695-8212 
htbanks@ncsu.edu 


20090319164 


AFRL-SR-AR-TR-09-0022 


The  public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information,  including 
suggestions  for  reducing  the  burden,  to  the  Department  of  Defense,  Executive  Service  Directorate  (0704-0188).  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law.  no 
person  shall  be  subject  to  any  penalty  for  failing  to  comply  with  a  collection  of  information  if  it  does  not  display  a  currently  valid  0M8  control  number. 

PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ORGANIZATION. 


14.  ABSTRACT 

In  this  10  month  funded  effort,  we  report  progress  on  the  development  of  methods  in  a  number  of  specific  areas:  stochastic  and  deterministic 
models  for  complex  networks  and  development  of  inverse  problem  methodologies  (generalized  sensitivity  functions). T  hese  efforts  are  part  of  our 
continuing  fundamental  research  program  in  a  modeling,  estimation  and  control  methodology( theoretical,  statistical  and  computational)  for  systems 
in  the  presence  of  major  model  and  observation  uncertainties. 


15.  SUBJECT  TERMS 

stochastic  and  deterministic  models  for  complex  networks;  inverse  problem  methodologies;  generalized  sensitivity  functions 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 
ABSTRACT 

18.  NUMBER 
OF 

PAGES 

19a.  NAME  OF  RESPONSIBLE  PERSON 

a.  REPORT 

b.  ABSTRACT 

c.  THIS  PAGE 

H.  T.  Banks 

U 

U 

U 

UU 

27 

19b.  TELEPHONE  NUMBER  (Include  area  code) 

919-515-3968 

Standard  Form  298  (Rev.  8/98) 

Prescribed  by  ANSI  Std.  Z39.18 
Adobt  Professional  7.0 


ABSTRACT 


In  this  10  month  funded  effort,  we  report  progress  on  the  development  of  methods  in 
a  number  of  specific  areas:  stochastic  and  deterministic  models  for  complex  networks  and 
development  of  inverse  problem  methodologies  (generalized  sensitivity  functions).  These  ef¬ 
forts  are  part  of  our  continuing  fundamental  research  program  in  a  modeling,  estimation  and 
control  methodology  (theoretical,  statistical  and  computational)  for  systems  in  the  presence 
of  major  model  and  observation  uncertainties. 


1  Brief  Overview  Summary 

Motivated  by  examples  in  viscoelasticity  and  biomedicine  such  as  in  [P5,P6,P7],  we  consider 
in  [P2]  the  problem  of  estimating  a  modeling  parameter  0  using  a  least  squares  criterion 
Jd(y ,  0)  =  (y(U)  —  f(U,  $))2-  We  take  an  optimal  design  perspective  in  [P2].  Our  general 

premise  (illustrated  via  examples  in  [P2] )  is  that  in  any  data  collected,  the  information 
content  with  respect  to  estimating  6  may  vary  considerably  from  one  time  measurement  to 
another,  and  in  this  regard  some  measurements  may  be  much  more  informative  than  others. 
We  propose  mathematical  tools  which  can  be  used  to  collect  data  in  an  almost  optimal  way, 
by  specifying  the  duration  and  distribution  of  time  sampling  in  the  measurements  to  be 
taken,  consequently  improving  the  accuracy  (i.e.,  reduce  the  uncertainty  in  estimates)  of  the 
parameters  to  be  estimated. 

We  recall  the  concepts  of  traditional  and  generalized  sensitivity  functions  and  use  these 
to  develop  a  strategy  to  determine  the  “optimal”  duration  time  T  for  an  experiment;  this 
is  based  on  the  time  evolution  of  the  sensitivity  functions  and  of  the  condition  number  of 
the  Fisher  information  matrix.  We  illustrate  the  role  of  the  sensitivity  functions  as  tools 
in  optimal  design  of  experiments,  in  particular  in  finding  “best”  sampling  distributions. 
We  also  present  a  new  optimal  design  criterion  which  is  based  on  the  idea  of  finding  the 
time  distribution  which  gives  the  best  discrete  approximation  for  the  integral  version  of  the 
Fisher  information  matrix.  Numerical  examples  are  presented  throughout  [P2]  to  motivate 
and  illustrate  the  ideas.  Theoretical  foundations  are  presented  in  an  appendix  of  [P2]. 

Alternative  but  related  methods  for  optimal  design  of  experiments  are  presented  and 
explored  in  [P4]. 

Stochastic  differential  equation  (SDE)  models  offer  one  formulation  for  introducing  un¬ 
certainty  in  human  interactions  in  a  dynamic  social  network  model  based  on  static  and/or 
deterministic  ordinary  differential  equation  (ODE)  models.  A  coupled  SDE  system  for  agent, 
characteristics  and  connectivities  was  developed  and  investigated  in  [7].  This  SDE  model 
(which  tacitly  assumed  instantaneous  influence  between  agents  with  connectivity)  may  be 
improved  by  including  delays  in  an  SDE  model  or  in  an  equivalent  Fokker-Planck  (FP) 
model  if  such  exists.  The  coupled  model  of  [7]  involved  discontinuities  and  did  not  yield  a 
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Markov  diffusion  process  (for  which  an  equivalent  Fokker-Planck  formulation  is  possible). 
In  this  effort  we  formulate  a  new  smooth  vector  SDE  system  and  demonstrate  that  it  gen¬ 
erates  a  Markov  diffusion  process  and  provides  computational  results  equivalent  to  those  of 
the  earlier  model  of  [7].  The  new  model  is  compared  to  the  earlier  model  via  simulations. 
We  derive  an  equivalent  Fokker-Planck  formulation  to  this  new  SDE  system.  Numerical 
methods  to  implement  the  FP  model  are  being  formulated  and  tested.  Versions  of  these 
equivalent  models  will  be  modified  first  to  include  deterministic  delays  and  then  delays  with 
some  stochasticiity / uncertainty. 

1.1  Publications  supported  in  part  under  this  grant 

[PI]  H.  T.  Banks  and  M.  Pedersen,  Well-posedness  of  inverse  problems  for  systems  with 
time  dependent  parameters,  CRSC-TR08-10,  August,  2008;  Arabian  Journal  of  Science  and 
Engineering:  Mathematics ,  to  appear. 

[P2]  H.T.  Banks,  Sava  Dediu,  Stacey  L.  Ernstberger  and  Franz  Kappel,  A  new  approach 
to  optimal  design  problems,  CRSC-TR08-12,  September,  2008;  Inverse  Problems ,  submitted. 

[P3]  H.T.  Banks,  J.E.  Banks  and  S.  L.  Joyner,  Estimation  in  time-delay  modeling  of 
insecticide-induced  mortality,  CRSC-TR08-15,  October,  2008;  J.  Inverse  and  Ill-posed  Prob¬ 
lems,  submitted. 

[P4]  H.T.  Banks.  J.L.  Davis,  S.L.  Ernstberger,  S.  Hu,  E.  Artimovich,  and  A.K.  Dliar, 
Experimental  design  and  estimation  of  growth  rate  distributions  in  size-structured  shrimp 
populations,  CRSC-TR08-20,  November,  2008;  Inverse  Problems ,  submitted. 

[P5]  H.T.  Banks  and  J.R.  Samuels,  Jr.,  Detection  of  cardiac  occlusions  using  viscoelastic 
wave  propagation,  CRSC-TR08-23,  December,  2008;  Advances  in  Applied  Mathematics  and 
Mechanics ,  submitted. 

[P6]  D.  Valdez-Jasso,  H.T.  Banks,  M.A.  Haider,  D.  Bia,  Y.  Zocalo,  R.L.  Armentano, 
and  M.S.  Olufsen,  Viscoelastic  models  for  passive  arterial  wall  dynamics,  CRSC-TR08-24, 
December,  2008;  Advances  in  Applied  Mathematics  and  Mechanics ,  submitted. 

[P7  H.T.  Banks,  Kathleen  Holm,  Nathan  C.  Wanner,  Ariel  Cintron- Arias,  Grace  M. 
Kepler  and  .Jeffrey  D.  Wetherington,  A  mathematical  model  for  the  first-pass  dynamics  of 
antibiotics  acting  on  the  cardiovascular  system,  CRSC-TR08-25,  December,  2008;  Mathe¬ 
matical  and  Computer  Modeling ,  submitted. 

[P8]  H.  T.  Banks.  K.  L.  Rehin  and  Karyn  Sutton,  Dynamic  social  network  models  with 
stochasticity  and  delays,  in  preparation. 

2  Personnel  Supported  (at  least  in  part  with  this  grant) 

•  H.T.  Banks,  Prof.,  North  Carolina  State  University 

•  G.  M.  Kepler,  Res.  Assoc.,  North  Carolina  State  University 

•  S.  Dediu,  Post  doc,  North  Carolina  State  University 
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•  S.  Hu,  Post  doc,  North  Carolina  State  University 

•  K.  Sutton,  Post  doc,  North  Carolina  State  University, 

•  S.  Joyner,  Graduate  Student  at  North  Carolina  State  University 

•  J.  Samuels,  Graduate  Student  at  North  Carolina  State  University 


3  Relevance  of  Research  to  DOD/AF  Missions 

The  development  of  general  classes  of  nonlinear  complex  nodal  network  models  with  in¬ 
herent  uncertainties  is  directly  relevant  to  several  recently  announced  “Discovery  Challenge 
Thrusts”  in  AFOSR-BAA-2007-08.  In  particular,  our  efforts  will  address  basic  research  is¬ 
sues  that  will  be  related  to  the  programs  announced  on  “  Robust  Decision  Making”  and 
“Complex  Networked  Systems”  which  are  concerned  with  operational  environments  where 
“unexpected  events,  attacks  on  and  degradations  of  the  network  centric  system  and  sud¬ 
den  changes  in  goals  and  objectives  ”  must  be  accommodated.  The  need  for  new  modeling 
paradigms  which  include  robustness  in  estimation  and  control  in  models  with  uncertainty  is 
well  recognized. 

We  expect  the  developed  methodologies  to  be  of  direct  interest  to  scientists  and  engineers 
in  the  Information  Directorate  at  AFRL,  Rome,  NY,  in  their  ongoing  efforts  on  information 
nodal  network  systems,  security  and  information  warfare. 


4  Detailed  Research  Summaries 

4.1  Generalized  Sensitivity  and  Optimal  Design 

We  consider  first  a  parameter  estimation  problem  for  the  general  nonlinear  dynamical  syst  em 

x{t)  =  g{t,x(t),9) 
x(t0)  =  x0, 


with  discrete  time  observations 


Vj  =  f(tj,  0)  +  tj  =  Cx(tj,  0)  +  tj,  j  =  l,...,n  (4.2) 

where  x,  g  €  RiV,  /,  tj  €  RAf  and  0  €  Rp.  The  matrix  C  is  an  MxJV  matrix  which  gives  the 
observation  data  in  terms  of  the  components  of  the  state  variable  x. 

Generalized  sensitivity  functions  (GSFs)  were  introduced  recently  by  Thoniaseth  and  Co- 
belli  [26  as  a  new  tool  in  identification  studies  to  analyze  the  distribution  of  the  information 
content  (with  respect  to  the  model  parameters)  of  the  output  variables  of  a  system  for  a 
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given  set  of  observations.  More  precisely,  the  generalized  sensitivity  functions  describe  the 
sensitivity  of  the  parameter  estimates  obtained  from  an  inverse  problem  with  respect  to  the 
data  measurements  used  in  the  inverse  problem. 

For  a  scalar  observation  model  with  discrete  time  measurements  (i.e.,  when  M  =  1  and 
C  is  a  1  x  N  array  in  (4.2)),  the  generalized  sensitivity  functions  (GSFs)  are  defined  as 


gs(*z)  =  E  ^TTT^1  x  W/(*i, 0o)]  •  Vtf/te. 0„), 

t=i  a 


(4.3) 


where  {£;},  l  =  1 ,n,  is  the  time  distribution  where  the  measurements  are  taken, 


F  =  E  0o)T  (4.4) 

is  the  corresponding  Fisher  information  matrix,  and  9q  is  the  usual  “true”  parameter  value 
tacitly  assumed  to  exist  in  standard  statistical  theories.  The  symbol  represents  element- 
by-element  vector  multiplication.  For  the  motivation  and  details  which  led  to  the  definition 
above,  the  reader  should  consult  [26]  and  [11].  The  Fisher  information  matrix  measures  the 
information  content  of  the  data  corresponding  to  the  model  parameters.  In  (4.3)  we  see  that 
this  information  is  transferred  to  the  GSFs,  making  them  appropriate  tools  to  indicate  the 
relevance  of  the  measurements  to  the  estimates  obtained  in  parameter  estimation  problems. 

We  note  that  the  generalized  sensitivity  functions  (4.3)  are  vector-valued  functions  having 
the  same  dimension  as  the  parameter  9  being  estimated.  The  fc-th  component  of  the  vector 
function  gs  represents  the  generalized  sensitivity  function  with  respect  to  9^.  The  GSFs 
are  defined  only  at  the  discrete  time  points  {tj,j  =  l,...,n}  and  they  are  cumulative 
functions,  at  each  time  point  ti  taking  into  account  only  the  contributions  to  estimates  of 
the  measurements  up  to  ti.  thus  indicating  the  influence  of  measurements  up  to  that  time 
on  the  parameter  estimates  obtained. 

From  (4.3)  it  is  readily  observed  that  all  the  components  of  gs  are  one  at  the  end  of 
the  total  time  interval,  i.e.,  gs (t^j)  =  1.  If  one  defines  gs(t)  =  0  for  t  <  t\  (gs  is  zero 
when  no  measurement  is  collected),  then  each  component  of  gs  varies  from  0  to  1  during 
the  experiment.  As  developed  in  26],  the  time  subinterval  during  which  this  transition 
has  the  sharpest  increase  corresponds  to  measurements  which  provide  the  most  information 
on  possible  variations  in  the  corresponding  estimated  model  parameters.  The  amount  of 
information  with  respect  to  a  parameter  9 *  is  directly  related  to  the  rate  of  change  of  the 
corresponding  GSF;  thus  sharp  increases  in  gs^,  indicate  a  high  rate  in  additional  information 
about  9^  being  provided  by  the  new  measurements  in  that  time  period. 

For  finite  dimensional  systems  ( N  <  oo)  with  a  finite  number  (p  <  oo)  of  parameters, 
the  numerical  implementation  of  the  generalized  sensitivity  functions  (4.3)  is  straightforward, 
since  the  gradient  of  /  with  respect  to  9  (or  rr0)  is  simply  the  Jacobian  of  x  with  respect  to 
9  (or  xo)  multiplied  by  C.  These  Jacobian  matrices  can  be  obtained  by  numerically  solving 
sensitivity  ODE  systems  coupled  with  the  system  (4.1). 
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We  have  effectively  used  GSFs  in  several  problems  [10,  4]  with  finite  dimensional  pa¬ 
rameters  in  ordinary  differential  equation  systems.  The  primary  contribution  of  generalized 
sensitivity  functions  is  that  they  can  indicate  regions  of  high  information  content  where,  if 
additional  data  points  are  taken,  one  can  generally  improve  the  existing  parameter  estimates. 
This  has  rather  obvious  applications  to  experimental  design.  Moreover,  by  visually  investi¬ 
gating  the  dynamics  of  GSF  curves,  one  can  potentially  identify  subsets  of  parameters  which 
are  highly  correlated.  While  there  are  still  some  heuristics  underlying  the  development  of 
GSFs  (the  presentation  in  the  appendiees  of  [5,  11]  offers  the  most  rigorous  presentation  to 
date),  it  is  clear  that  some  version  of  the  GSF  theory  should  become  a  valuable  tool  for  sci¬ 
entific  and  engineering  investigators  needing  to  estimate  parameters  in  dynamical  systems. 
Although  the  GSF  theory  provides  useful  information  in  parameter  estimation,  the  most 
insight  can  be  gained  when  the  GSFs  are  used  in  conjunction  with  traditional  sensitivity 
functions. 

In  order  to  illustrate  our  ideas,  we  changed  the  traditional  order  of  the  steps,  and  we 
assume  that  we  need  to  estimate  the  parameter  8  from  a  set  of  observations  of  the  output 
variable  y  =  Cx  which  are  not  known  a  prion .  Indeed,  let  us  assume  that  the  only  things  we 
know  for  the  moment  are  the  physical  problem  and  the  mathematical  model  (4.1)-(4.2)  used 
to  model  it,  but  we  have  the  freedom  to  instruct  the  experimentalists  as  to  the  number  of 
observations  to  take,  and  more  importantly  at  what  times  to  take  them.  We  also  assume  that 
the  underlying  model  (4.1)-(4.2)  is  a  good  representation  of  the  physical  problem  (i.e.,  there 
is  agreement  between  the  observed  data  and  the  corresponding  values  generated  by  (4.1) 
(4.2)  with  an  appropriate  choice  of  0),  and  that  we  have  a  criterion  to  ascertain  the  quality 
of  the  estimates  obtained  (for  instance,  the  standard  errors).  Given  these  assumptions,  it 
makes  sense,  once  we  have  a  mathematical  model  and  an  algorithm  to  perform  the  inversion 
(least  squares  for  example),  to  look  for  a  strategy  for  collecting  experimental  data  in  the 
best  possible  way.  We  want  to  take  maximum  advantage  of  the  mathematical  structure  of 
the  model  in  order  to  indicate  the  optimal  number  of  observations  to  take,  as  well  as  their 
optimal  temporal  distribution,  thus  obtaining  better  accuracy  for  the  estimates. 

Some  of  the  most  important  questions  in  parameter  estimation  problems  when  considering 
the  dynamical  system  (4.1)-(4.2)  are: 

1.  How  do  we  choose  the  duration  T  of  the  experiment,  or  the  interval  [0 ,T]  from  where 
to  sample  data,  in  order  to  obtain  the  best  estimates  for  0 ? 

2.  Once  the  duration  T  is  established,  what  is  the  optimal  sampling  distribution  of  an  a 
priori  given  number  of  measurements  in  the  interval  [0,  T],  in  order  to  obtain  the  most 
accurate  estimates? 

The  first  question  hints  to  the  minimum  time  interval  which  contains  most  of  the  relevant 
information  with  respect  to  the  estimation  of  8.  The  second  question  constitutes  the  main 
topic  of  the  optimal  experimental  design  literature,  and  the  traditional  way  to  answer  it  is 
to  find  time  distributions  which  optimize  some  specific  “design”  functions.  In  general,  an 
optimal  design  minimizes  or  maximizes  a  given  function  of  the  Fisher  information  matrix. 
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However,  this  matrix  depends  on  the  true  parameter  80 ,  which  in  practice  is  unknown.  In 
order  to  overcome  this  difficulty  one  usually  substitutes  an  initial  guess  9*  for  The 
resulting  designs  are  called  locally  optimal  designs.  Of  the  numerous  design  strategies  found 
in  the  literature  (for  example,  see  [12,  17,  18,  24])  we  briefly  enumerate  here  the  most  popular, 
which  are: 

a)  the  D-optirnal  design ,  which  arguably  maximizes  the  determinant  det F(0*)  of  the 
Fisher  information  matrix.  It  is  largely  accepted  in  the  literature  because  of  its  ap¬ 
pealing  geometrical  interpretation:  the  asymptotic  confidence  regions  for  the  maximum 
likelihood  estimate  of  6  are  ellipsoids.  One  can  show  that  an  optimal  design  with  re¬ 
spect  to  this  criterion  yields  approximately  a  minimal  value  for  the  volume  of  the 
tolerance  ellipsoid  of  the  estimates; 

b)  the  c- optimal  design ,  which  minimizes  the  variance  of  a  linear  combination  g(9)  =  cT6 
of  the  parameters  to  be  estimated.  Up  to  a  multiplicative  constant,  the  asymptotic 
covariance  of  the  maximum  likelihood  of  g  is  given  by  cT F~l(9*)c,  and  an  optimal 
design  with  respect  to  this  criterion  minimizes  this  variance.  In  the  particular  case 
when  c  is  the  z-tli  unit  vector,  i.e.,  c  =  e*,  the  c-optimal  design  minimizes  the  variance 
of  the  least  squares  estimate  for  the  z-tli  parameter  Of, 

c)  the  E -optimal  design ,  which  maximizes  the  smallest  eigenvalue  Am{n(F(0*))  of  the 
Fisher  information  matrix.  This  is  equivalent  to  minimizing  the  maximum  eigenvalue 
Ama X(F  1(9*))  of  the  inverse  F~l  (the  covariance  matrix),  or  to  minimizing  the  worst 
variance  among  all  estimates  c*T 8  for  the  linear  combinations  cJ8  with  a  given  norm 
cTc  —  1.  Geometrically,  the  jF-optiinal  design  minimizes  the  maximum  diameter  of  the 
asymptotic  confidence  ellipsoids  for  9. 

The  practical  importance  of  the  issues  addressed  by  the  previous  two  questions  is  obvious 
for  those  experiments  where  the  cost  of  measurements  (invasive  procedures,  expensive  tech¬ 
nology,  necessity  of  highly  qualified  personnel,  etc.)  is  high  and  where  one  wants  to  avoid 
running  the  experiments  longer  than  necessary.  These  questions  suggest  that  we  investigate 
the  temporal  distribution  of  the  information  content  with  respect  to  parameters  we  want  to 
estimate,  the  sensitivity  of  the  model  output  and  of  the  parameter  estimates  with  respect  to 
model  parameters,  and  lead  us  to  introduce  mathematical  tools  which  to  describe  them. 

Our  strategy  in  answering  the  previous  questions  is  to  assume  that  we  have  continuous 
time  measurements  in  [0,  T)  and  develop  tools  to  quantify  the  information  content  with 
respect  to  9  (imbedded  in  rf)  throughout  the  interval  [0 ,  T].  In  order  to  carry  this  out,  wc 
reformulate  the  estimation  problem  by  using  an  integral  form  of  the  least  squares,  integral 
optimality  conditions,  and  integral  version  of  the  Fisher  information  matrix  (its  integral  form 
on  the  interval  [0,  T]).  First  wc  show  that  the  traditional  and  generalized  sensitivity  functions 
(GSF)  can  indicate  regions  of  high  information  content,  where  if  additional  data  points  are 
sampled,  the  corresponding  parameter  estimates  are  improved  (as  wc  have  noted  GSF  were 
introduced  first,  to  our  knowledge,  by  Thomaseth  and  Cobelli,  see  [4,  11,  26]).  Therefore 
we  advocate  their  utilization  as  new  tools  in  experimental  design.  In  [5]  wc  propose  a  new 
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optimal  design  strategy  which  is  based  on  the  idea  of  finding  a  discrete  sampling  distribution 
on  [0,  T],  which  best  approximates  the  integral  version  of  the  Fisher  information  matrix  with 
a  discretized  version.  To  carry  this  out,  we  formulate  in  [5]  a  new  continuous  GSF  based 
on  probability  distributions  representing  the  sampling  which  can  then  be  optimized  over  the 
metric  space  of  probability  measures  formulated  with  the  Prohorov  metric  [2,  23]. 

4.2  Nodal  Networks  with  Uncertainty 

Networks  with  uncertainty  are  ubiquitous  in  science  and  engineering,  but  especially  in  com¬ 
plex  nonlinear  systems  of  interest  to  DOD.  In  particular,  nodal  networks  with  uncertainty  are 
fundamental  to  the  understanding  of  problems  involving  dynamic  resource  management  and 
allocation.  Examples  range  from  service  scheduling  and  material  supply  chains  to  produc¬ 
tion  systems  and  modern  high  speed  communication  networks  where  information  is  routed 
between  nodes  in  a  distributed  information  system  with  a  global  grid. 

Social  network  structures  can  be  represented  by  graphs  in  wdiich  individuals,  or  agents, 
are  represented  by  vertices  and  the  connections  between  agents  are  represented  by  edges. 
Agents  are  defined  by  measurable  characteristics ,  and  the  connections  between  agents  can  be 
defined  as  either  unidirectional,  thus  creating  a  directed  graph,  or  bidirectional,  thus  creating 
an  undirected  graph.  In  social  network  graphs,  the  relationship  between  two  agents,  or  a 
pairwise  relationship,  is  considered  to  be  the  most  basic  structure  used  in  constructing  the 
graph.  More  complicated  social  models,  even  triads  including  three  agents,  are  created  from 
using  multiple  pairwise  relationships.  Research  on  social  networks  may  be  applied  to  many 
areas,  including  social  group  structure,  information  sharing,  and  disease  outbreak. 

Interactions  between  agents  often  have  uncertain  and  unpredictable  results,  and  many 
static  and  deterministic  social  network  models  fail  to  consider  variability  and  delay.  Complex 
nonlinear  systems  are  more  suited  for  modeling  the  uncertainty  found  in  human  interactions. 
A  dynamic  social  network  model  that  includes  delay  and  accounts  for  variability  could  allow 
governmental  organizations  to  predict  the  propagation  of  information  through  a  terrorist  net¬ 
work  2],  allowr  local  officials  to  examine  how  a  communicable  disease  would  spread  through 
their  community  [20],  or  enable  corporations  to  simulate  scenarios  in  wdiich  their  supply 
chain  operates  under  strained  conditions  or  with  perturbations  [10]. 

One  area  of  current  research  in  social  networks  is  dynamics  of  connections  and  agents 
over  a  time.  The  model  first  proposed  by  [7]  utilizes  a  coupled  system  of  stochastic  dif¬ 
ferential  equations  (SDK’s).  Agents  are  defined  by  quantitatively  described  characteristics 
and  the  strength  of  the  pairwise  relationship,  or  degree  of  connectivity ,  joining  two  agents 
is  determined  by  these  characteristics.  It  is  assumed  that  agents  participate  in  hornophilic 
relationships  and  that  two  agents  are  more  likely  to  have  a  positive  degree  of  connectivity 
if  their  characteristics  are  similar.  The  proposed  model  accounts  for  changes  in  agents  and 
relationships  through  both  observable  interactions  and  unobservable  events  represented  as 
noise. 

Another  factor  in  social  network  dynamics  that  can  be  considered  is  delay.  The  effect 
one  agent  has  on  another  often  is  not  instantaneous;  like  diseases  in  a  body,  new  ideas  and 
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opinions  are  subject  to  a  propagation  time  in  an  agent’s  rnind  before  the  agent  accepts  or 
rejects  the  new  concept.  An  agent  does  not  pass  on  ideas,  objects,  and  diseases  immediately 
after  they  have  come  in  contact  with  or  possession  of  them.  Consider  a  conversation  between 
two  people  of  differing  religious  views.  It  is  unlikely  that  one  will  immediately  convert  to 
the  other’s  faith;  however,  after  repeated  contact  and  reflection  by  each  individual,  the  two 
may  become  more  understanding  of  each  others’  religions.  Adding  constant  and  stochastic 
delays  to  the  proposed  model  will  result  in  a  more  realistic  simulation  of  agent  interactions 
and  internal  changes  in  the  agents. 

The  proposed  model  is  built  from  the  model  created  by  [7].  The  model  is  used  with  the 
assumption  that  time,  degrees  of  connectivity,  and  agent  characteristics  are  continuous.  The 
coupled  SDE  model  proposed  in  [7]  was  revisited  and  used  to  reproduce  previous  results. 
The  SDE’s  were  then  changed  so  that  the  coefficients  used  in  the  model  are  smooth,  allowing 
the  model  to  satisfy  the  conditions  of  a  diffusion  process,  and  the  results  of  the  model  with 
smooth  coefficients  are  compared  to  the  results  of  the  original  model  to  verify  that  the  change 
in  coefficients  does  not  affect  the  visible  behavior  of  the  model.  These  equations  arc  then 
converted  to  an  equivalent  Fokker-Planck  (FP)  model.  Numerical  implementation  of  the  FP 
model  as  an  alternative  modeling  and  numerical  approach,  with  attention  to  ability  to  add 
delays  to  the  system,  is  under  investigation. 

4.2.1  Basic  Stochastic  Differential  Equation  Model 

Denote  the  vector  of  characteristics  of  agent  z,  z  £  M  =  1, . . . ,  N  at  time  t  by  C*(t)  £ 
K  C  Rm,  where  K  is  a  compact  constraint  set  for  the  values  of  characteristics  and  m  is  the 
number  of  characteristics.  Refer  to  the  degree  of  connectivity  between  agent  z  and  agent  z' 
at  time  t  by  e(z,z',£)  €  R  such  that  an  agent’s  connectivity  to  itself  e(z,z,£)  =  0  Vi,  t  and  an 
edge  between  z  and  if  exists  if  and  only  if  e(z,  z',  t)  >  0.  The  strength  of  the  link  between 
agent  z  and  agent  z'  is  not  bidirectional,  so  it  is  possible  that  e(z,z',£)  ^  e(z',z,t).  Define 
Ai(t)  =  {z'  £  J\f  :  e(z,  z',  t )  >  0}  to  be  the  set  of  all  agents  to  which  agent  i  is  linked  at  time 
f,  and  let  \Ai(t)\  be  the  number  of  elements  in  Ai(t). 

The  system  of  coupled  stochastic  differential  equations  as  defined  by  [7]  used  as  the 
starting  point  in  this  model  are 


dCi(t) 


de(i,  i',  t ) 


A 


AAt)  I 

n  '  1  i'eAi(t) 


X;  fc,(t)  -  Ci(t)]dt  + 


(4.5) 


(4.6) 


In  (4.5),  the  constant  j3{  is  an  agent-specific  value  that  represents  the  desire  of  agent  z 
to  be  more  like  agents  to  which  it  is  linked  and  crWp(-)  is  a  Wiener  process  with  variance 
a2.  In  (4.6),  /(£)  =  2e“^  —  1  where  the  constant  b  >  0  controls  the  rate  at  which  edges  arc 
formed  and  ||  •  ||  denotes  the  square  of  the  vector  norm  in  Mm,  and  7dW^,(*)  is  a  Wiener 
process  with  variance  72. 
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4.2.2  Stochastic  Differential  Equation  with  Smooth  Coefficients 

The  model  described  by  (4.5)  does  not  fit  the  definition  of  a  usual  diffusion  process  because 
the  values  of  |  Al(t)  |  change  discontinuously.  An  alternative  model  is  proposed  in  [7]: 


dCl(t)  =  E ^fe}ge(i’A,) 


C i'(t)  -  C i(t)  dt  +  ocTWf(t) 


(1  7) 


The  sum  e(z,  i\  t)  is  defined  for  all  sums  of  all  possibles  values  of  e(? ,  z',  t )  and  thus 
is  continuous.  While  it  has  smooth  coefficients,  this  model  fails  to  follow  the  assumption  of 
homophily  and  enables  agents  who  have  a  non-positive  connectivity  to  a  particular  agent 
to  influence  the  characteristics  of  that  particular  agent.  This  in  turn  causes  the  degrees  of 
connectivity  e(z,  i',  t(n))  — >  oo  as  n  — >  oo,  resulting  in  single  cluster  scenarios  for  any  value 
of  b  in  the  deterministic  case,  as  illustrated  by  Figure  1. 


Degree  of  connectivity  e(IO.i'.t)  for  b=0.02 


Figure  1:  Degree  of  connectivity  e(10,  z',  i),  z7  ^  10  of  agent  10  with  other  agents  for  b  =  0.02 
using  smooth  model  proposed  in  [7]. 

Consider  the  agents  i\  i 9  ^  2,  7  and  agent  2  in  the  case  study  used  in  [7]  in  which 
C2/(0)  =  (10  10)  and  C2(0)  =  (  —  10,  —10).  Assume  that  these  agents’  characteristics  remain 
the  same  for  several  time  steps  in  the  simulation  while  the  degrees  of  connectivity  e(?’V2,£) 
and  e(2,z',£)  become  negative  numbers  so  that  e(2,  i\tn)  =  —1.5  at  some  time  tn.  Ignoring 
agent  7,  (4.7)  would  reduce  in  the  deterministic  case  to 
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dC2(tn)  = 

dC2(tn)  = 
dC2(tn)  = 


= - S-— -  Y,  «0.  «.)  fc,.((„)  -  C2(U 

2^'/2  e\Zi  1  ’  L 

+<riWp(i„) 

-52-1.5 

-19^ 


eft 


12 

20 

20 


i'^2 


dt 


no' 

1 - 

0 

T— 1 

1 

1 _ 

[io_ 

— 10J 

dt  +  0dW?(tn) 


As  noted  before,  agents  to  which  agent  2  has  a  negative  degree  of  connectivity  should  not 
alter  the  characteristic  values  of  agent  2,  but  this  rule  is  clearly  violated  by  this  proposed 
smooth  model. 

To  avoid  agents  affecting  each  other  at  values  of  e(i,i\t)  <  0,  several  augmentations  to 
the  sum  c(i,  t)  were  considered.  One  proposed  change  that  enabled  the  model  to 

produce  results  similar  to  those  produced  by  (4.5)  was 


dC.it)  - 


Pi 


^  *v*  ®(^i 


i\t) 


E 

•Vi 

•'€^i(f) 


e(z,z",i)  C*(0-Ci(*)  di  +  arfWp(«)  (1.8) 


Like  in  (4.5),  At(£)  is  the  set  of  agents  i!  /  i  that  have  connectivity  >  0.  As 

the  creation  of  Ai{t)  requires  additional  calculations  arid  computer  resources,  this  method  is 
less  efficient  than  the  ideal  method.  Imposing  the  additional  condition  on  the  sum  may  also 
have  unwanted  effects  when  this  equation  is  converted  to  a  Fokkcr-Planck  equation.  The 
SDK  system  with  smooth  coefficients  containing  an  equivalent  method  of  expressing  (4.8)  is 


dC.it) 


rie(z,  t) 


|Q/(t)  -  Cjjt)  dt  +  adWfit) 
/(||Ci-Cl,||2)eft  +  7rfWJt<,(«). 


(4.9) 

(4.10) 


In  (4.8),  p  =  -  (e(z,  i\  t)  T  |e(z,  i\  t)|)  As  no  agents  i*  are  to  affect  agent  i  if  e(L  /•',  t)  <  0. 
it  is  assumed  that  ^  —  -  =  0  if  JT,  =  0.  Thus  this  model  will  not  fail  if  single-agent 

clusters  arise  during  a  simulation.  This  model  maintains  the  requirement  of  homophily 
between  agents  and  prevents  agents  with  11011-positive  degrees  of  connectivity  affecting  other 
agents.  When  used  to  run  simulations  of  the  test  cases  first  run  with  (4.5),  results  similar 
to  those  of  the  earlier  model  are  produced. 


4.2.3  Numerical  Scheme  for  Stochastic  Models 

We  continued  to  use  the  stochastic  analogue  of  the  classical  fourth-order  Runge-Kutta 
method  to  solve  the  system  of  differential  equations  as  implemented  by  [7].  Approximating 
dCiit )  remains  nearly  the  same: 
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C  i(tn)  —  Ci(tn-i)  +  poFoh  +  p\Fih  +  P2F2I1  +  P3F3I1  +  a  AW^( 


(4.11) 


where 


i\^n—  1 


Cta)(4-i)  =  +  Ifo/H- 1(tAW® 


Fi  —  g  (tn-  1  +  ^11 


C?\tn-X) 

f2 


Q(tn-i)  +  2^*^  +  — crAW^ 


F,  = 


g(*n-l  +  l/l.Cf^n-i 

Ci(^n-l)  +  F2/1  +  (TAWjJ 

g^n-1  + 


Po  =  P3  =  7:  and  Pi  =  P2  =  7:. 

6  3 

The  value  of  <7  was  added  to  the  calculation  of  Cz-^(tn_i),  Cp^(tn_i),  and  Cp^(tn_i) 
to  scale  the  value  of  AW^  so  that  it  reflected  the  value  of  aAW^  used  in  (4.11).  The 
definition  of  g(t,  C i(t))  is  depedent  upon  the  model.  For  approximation  of  the  basic  SDE  , 

g (t,  Ci(t))  =  t-Avt  £\,€i4  m  Ci/(t)-C<(t)  ,  and  for  the  smooth  coefficient  SDE,  g(t,  C ,(£))  = 


Equation  (4.6)  may  be  solved  using  a  simplified  manner  since  its  result  is  not  dependent 
upon  the  value  of  the  degree  of  connectivity  e(z,  i',i).  This  can  be  expressed  as 

e(i,i',tn)  =  e(i,i',tn)  +p0/(||Q(t„_1)  —  C</(*n_i)||2)/i 

+(pi  +P2)/(||Cj(fn_i  +  -fl)  —  Cj' (tn— 1  +  2OH 

+P3/(||Q(tn)  -Cay||2)/*  +  7AW^.  (4.12) 

Here  AW-n  —  W^,(tn)  —  W eiif(tn-\)  are  standard  Wiener  increments  with  the  distribu¬ 
tion  N(0,  h)  where  h  =  A tn.  We  assume  that  the  value  of  C i(tn)  +  ^ h  can  be  approximated 
linearly  from  Q(in_i)  to  C i(tn)  when  A tn  is  sufficiently  small  and  therefore  C +  -/i) 
may  be  approximated  by  cd^-i)+c1(^ri) 
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4.2.4  Numerical  Simulations 


We  use  the  same  case  study  concerning  10  agents  as  in  [7]  over  a  time  interval  of  t  G  [0, 10]. 
To  avoid  being  undefined,  we  calculate  e(i, i',0)  =  /(J|Q(0)  —  Q/(0)||2^/i.  This  will 

cause  all  connectivities  to  start  with  a  nonzero  value  that  is  also  indicative  of  the  sign  of  the 
degree  of  connectivity  between  two  agents  as  determined  by  their  characteristics. 

4.2.5  Deterministic  Case 

The  reproduction  of  the  existing  model  performs  similarly  to  that  of  the  paper  [7]  in  the 
deterministic  case  (when  a  =  0  and  8  =  0).  By  properly  implementing  (4.6),  the  values  of 
b  that  determined  the  number  of  clusters  created  after  the  same  time  as  used  in  [7]  were 
found  to  be  different.  One  cluster  formed  when  b  G  (0,0.00155);  two  clusters  formed  when 
b  G  [0.00155,0.01390);  and  three  clusters  formed  when  b  G  [0.01390,  oo).  Two  agents  are 
members  of  different  clusters  if  their  connectivities  <  0  and  e(i',i,  t)  <  0.  The 

clusters  are  identical  to  those  described  in  [7]:  the  one  cluster  scenarios  had  a  result  of  one 
cluster  of  all  10  agents;  the  two  cluster  scenarios  had  a  result  of  one  cluster  of  agents  2  and 
7  and  a  second  cluster  of  all  other  agents;  and  the  three  cluster  scenario  had  a  result  of  one 
cluster  of  only  agent  2,  one  cluster  of  only  agent  7,  and  one  cluster  of  all  other  agents,  and 
the  degrees  of  connectivity  between  agents  changed  at  a  rate  identical  to  the  rates  found  in 

[7]- 


Degree  of  connectivity  e(IO.i’.t)  for  b=0  02 


Figure  2:  Degrees  of  connectivity  e(10,  i' ,  t),i'  7^  10  of  agent  10  with  other  agents  for  b  —  0.02 
(three  cluster  scenario). 
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Degree  of  connectivity  e(10,i’,t)  for  b=0.005 


Figure  3:  Degree  of  connectivity  e(10,  i\  i),  if  ^  10  of  agent  10  with  other  agents  for  b  =  0.005 
(two  cluster  scenario). 


Degree  of  connectivity  e(10,i',t)  for  b=0.000775 


Figure  4:  Degree  of  connectivity  e(10,  7^  10  of  agent  10  with  other  agents  for  b  = 

0.000775  (one  cluster  scenario). 
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Noise  Enriched  Situations 


Figure  5:  The  (<r,  7)  noise-enriched  surface  along  a  £  [0,40],  7  £  [0, 12]  plotted  as  a  number 
out  of  100  simulations. 

4.2.6  Stochastic  Case 

Four  regimes  are  defined  in  [7]  into  which  the  results  of  a  simulated  trial  may  be  organized. 
These  four  regimes  are  labeled  essentially  deterministic,  noise  enriched,  noise  enlarged  (with 
two-  and  three-cluster  possibilities),  and  noise  dominated. 

To  investigate  the  occurrence  of  the  regimes  defined  in  [7]  using  a  constant  b  and  variable 
a  and  7  we  ran  simulations,  with  time  step  h  =  Atn  =  0.05,  for  selected  values  of  a  £  [0, 40] 
and  7  £  [0, 12  with  fixed  b  =  0.000775.  We  carried  out  100  simulations  for  each  pair  (07..  7 /) 
in  a  partition  with  values  0*  —  k  =  0, 1, . . . ,  40;  7/  =  0.25/,/  =  0,1,...,  48.  To  further 
investigate  the  relationship  between  <7,  5,  and  the  number  of  trials  that  produced  essentially 
deterministic  results,  we  also  carried  out  100  simulations  for  each  pair  (0^,7 /)  in  a  partition 
with  values  =  k  =  0, 1, ... ,  40;  7 /  =  0.0025/,  /  =  0,1,...,  48.  The  trends  exhibited  in  [7] 
Figures  10  -  13  are  also  present  in  the  charts  produced  by  the  proposed  model  albeit  on  a 
larger  domain  of  a  and  7.  Again,  the  surface  of  the  noise-dominated  regime  is  almost  the 
complement  of  the  surface  of  the  noise-emiched  regime  with  the  exception  of  the  occurrence 
of  some  essentially  deterministic  cases  at  low  values  of  a  and  7.  This  indicates  that  it  is 
highly  unlikely  that  the  inclusion  of  noise  in  the  model  will  cause  the  simulation  to  reach  a 
two-  or  three-clustered  scenario  if  a  b  in  the  range  of  values  that  will  produce  a  one-clustered 
scenario  is  selected. 
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Noise  Dominated  Situations 


Figure  6:  The  (cr,  7)  noise-dominated  surface  along  a  G  [0, 40],  7  G  [0, 12]  plotted  as  a  number 
out  of  100  simulations. 


Deterministic  Situations 
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Figure  7:  The  (cr,  7)  essentially  deterministic  surface  along  a  G  [0,40],  7  G  [0,0.12]  plotted 
as  a  number  out  of  100  simulations. 
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Noise  Enriched  Situations 


Figure  8:  The  (a.  7)  noise  enriched  surface  along  a  €  [0, 40],  7  €  [0, 0.12]  plotted  as  a  number 
out  of  100  simulations. 
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4.2.7  Stochastic  Model  with  Smooth  Coefficients 

The  coupled  SDE  model  with  smooth  coefficients  were  also  used  to  simulate  the  case  denoted 
above  and  produced  visually  identical  results.  In  the  deterministic  case,  the  values  of  b  that 
control  clustering  behavior  are  identical  to  those  found  in  the  deterministic  case  of  the 
basic  SDE  system;  that  is,  one  cluster  formed  when  b  G  (0,0.00155);  two  clusters  formed 
when  b  G  [0.00155,0.01390);  and  three  clusters  formed  when  b  G  [0.01390,  oo).  The  sample 
distributions  of  regimes  also  appear  to  be  very  similar  to  those  established  in  [7].  Instead 
of  using  agent  10  to  represent  all  agents  j  with  Cj(0)  =  (10, 10),  agent  1  was  used  with  no 
impact  of  the  behavior  of  the  model. 

To  verify  that  the  model  exhibited  random  behavior  when  a  ^  0  and  7  7^  0,  a  simulation 
of  the  one-cluster  scenario  was  run  using  the  parameters  a  =  5  and  7  —  1.  The  general 
behavior  exhibited  by  the  deterministic  case  is  also  evident  in  this  stochastic  case.  Note  that 
the  trend  of  increasing  connectivities  in  a  one-cluster  deterministic  case,  pictured  in  Figure 
4.2.7,  is  also  reflected  in  the  one-cluster  stochastic  case  pictured  in  Figure  4.2.7.  Also,  the 
tendency  for  characteristic  values  in  the  one-cluster  deterministic  case  to  approach  a  certain 
value  between  all  the  characteristic  values,  shown  in  Figure  4.2.7,  is  also  present  in  the 
one-clustcr  stochastic  case,  as  illustrated  by  Figure  4.2.7. 


Degree  of  connectivity  e(1  ,i‘,t)  for  b=0  000775 


Figure  9:  Degree  of  connectivity  e(l,  i\  £),  i'  7^  1  of  agent  1  with  other  agents  for  b  —  0.000775, 
a  =  0,  and  7  0  using  the  smooth  coefficient  model  (one  cluster  scenario,  deterministic 

case). 
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Characteristics  of  agents  1,  2,  and  7  over  time 
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Figure  10:  Characteristic  values  of  agents  1,  2,  and  7  over  t  G  [0, 10]  for  b  =  0.000775,  a  =  0, 
ami  7  =  0  (one  cluster  scenario,  deterministic  case). 


Noise  Enriched  Situations 


Figure  11:  The  (it.  7)  noise  enriched  surface  along  a  G  [0, 40],  7  G  [0, 12]  plotted  as  a  number 
out  of  10  simulations  using  the  smooth  coefficient  model. 
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Noise  Dominated  Situations 


Figure  12:  The  (a,  7)  noise  dominated  surface  along  a  G  [0,40]  7  G  [0,12]  plotted  as  a 
number  out  of  10  simulations  using  the  smooth  coefficient  model. 


Degree  of  connectivity  e(1  ,i',t)  for  b=0  000775 


Figure  13:  Degree  of  connectivity  e(l,i/,4),i/  ^  1  of  agent  1  with  other  agents  for  one  trial 
using  the  parameters  b  =  0.000775,  a  =  5,  and  7  =  1  using  the  smooth  coefficient  model 
(one  cluster  scenario,  stochastic  case). 
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Figure  14:  Characteristic  values  of  agents  1,  2,  and  7  over  t  E  [0, 10]  for  one  trial  using  the 
parameters  b  =  0.000775,  a  =  5,  and  7=1  (one  cluster  scenario,  stochastic  case). 
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4.2.8  Fokker-Planck  Equivalent  Formulation 

With  the  establishment  of  the  smooth  SDE  model  and  its  faithfulness  to  the  SDE  model 
first  proposed  by  [7],  the  development  of  an  equivalent  FP  model  may  proceed.  Let  {A"(£)} 
be  a  stochastic  process,  with  t  6  [to,  oo),  and  X(t)  a  random  variable  continuous  in  time  and 
state,  X(t)  £  (  — oc,  op).  Let  X(t)  satisfy  the  Ito  SDE 

dX(t)  =  a(X(t),  t)dt  +  f3(X(t),  t)d\V(t),  (4.13) 

where  d\V(t)  represents  an  infinitesimal  Wiener  increment,  and  W(t)  is  the  standard  Wiener 
process.  The  mean  and  variance  of  X(t)  is  then  a(X(t),t)  and  ft2(X(t),  i),  respectively. 
Conditions  (4.14)  -  (4.16) 


lim 

0  At  J 

r  oo 

/  1  y- 

—oo 

x\sp(y ,  t  +  At;  x,  t)dy  =  0, 

6  >2 

(4.11) 

lim 

A£-+0  At  j 

roo 

/  (y- 

-oo 

*)p(y,  t  +  At;  x,  t)dy  =  a(x,  t), 

(4.15) 

lim  — 

A*— 0  At  J 

—  oo 

x)2p(y,  t  +  At;  x,  t)cly  =  f 32{x ,  t), 

(4.16) 

establish  the  smoothness  of  coefficients  a(X(t),t)  and  /3(X(t),t).  Should  these  conditions 
be  satisfied,  it  can  be  shown  via  the  application  of  the  Chapman-Kohnogorov  equation  and 
Taylor  expansion  of  a  test  function,  that  X(t)  is  a  diffusion  process.  Equivalently,  this  means 
that  X(t)  satisfies  the  Kolmogorov  equations  with  drift  coefficient  a  and  diffusion  coefficient 
f32.  The  corresponding  forward  Kolmogorov  (or  Fokker-Planck)  equation  is  then 

dp  d(a(x,t)p)  ld2((32(x,t)p) 

U= - d^  +  2 - d? - '  (417) 

In  equation  (4.17),  p  is  the  pdf  of  the  stochastic  process  X(t),  and  the  transition  from 
state  x  at  time  t  to  state  y  at  time  s  is  then  p(y,  s ;  x,  t). 

Generalizing  (4.17)  to  M  dimensions  as  in  [22],  we  find  the  forward  Kolmogorov/Fokker- 
Planck  is  given  by 


dp 

dt 


where  X  =  (X\, . . . ,  X a/)  is  a  stochastic  process  such  that 


and 


(4.18) 


S'j  =  ^toE[AX*®AXM\Xi(t)  ~  XiMt)  =  x ,}. 
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If  we  let  Cf(f)  represent  the  kth  element  of  C,(f),  and  apply  (4.18),  the  forward  Kol¬ 
mogorov  formulation  corresponding  to  model  (4.9)  -  (4.10)  is  then 


dp 

dt 
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1=1 
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dC } 
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E*  «(*»*',  0 
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E*  e(M'.0 
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TV  JV 
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(01  p 
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+  EE^Mic'W-c‘(*)ii2M 


i=i  <-i  9<v 

ifsVp)  fayp)  ,  A  3V2p)  ffW) 

'  2  ^  a(c/)2  +  £  d(c?)2  -  ^a(CD2  1  ■ 


Here  p  is  the  joint  transition  probability  density  function  of  p(C,(£),  e(z,  i',  Summing 

the  drift  of  C*  ( 1 )  and  the  variances  of  the  Wiener  processes  for  each  characteristic  results  in 


dp 

dt 


N  m 
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EE9Cfc 


i=l  k= 1 
jV  /v 
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+  E  E  [/(IIc.'<0  -  c<(i)ll2);'] 
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,  1  V—  <J2(.'T~j,>  V—  V-  d2(l2I>) 

2  tr  fe  8<c‘)!  ^  ^ 


ae«' 


(4.19) 


4.2.0  Finite  Difference  Scheme  for  Fokker-Planck  Model 

The  finite  difference  approximation  formulated  by  [14]  preserves  the  properties  of  the  FP 
model.  This  numerical  scheme  yields  a  non-negative  solution  for  all  values  of  x,  conserves 
the  area  of  the  p.d.f.  p,  and  is  known  to  be  relatively  stable.  Additionally,  this  method 
allows  mesh  points  to  be  more  widely  spaced  than  in  some  competing  methods.  Define  the 
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following  functions 


G(x,  t) 


F(x,t) 


N 


EEfiLE*''1)  [<$(«>  -  c?(*)] 

1=1  k=  1  V  1  ’  }  <t>i 

+  5Z5Z2  [exp(-fe||Ci(f)  -Ci/(f)||2)  -  1] 


i=l  i;=l 

cr 


9  N  m  q  9  N  N 

^tEE5^4SE 


dp 


2  2  £^de(M',0 


(4.20) 


(4.21) 


For  any  x  G  x,  let  Ax  =  (xmax  —  xrnin)/R,  where  xrnax  and  xm*n  are  the  upper  and  lower 
boundaries  of  the  domain  of  x  G  x  and  R  is  the  number  of  intervals  into  which  the  domain 
of  x  shall  be  divided.  Let  At  =  T/S ,  where  S  is  the  number  of  time  steps.  The  mesh  points 

shall  be  given  by  xr  —  xrnin  +  r  Ax,  r  =  0, 1, 2, - /?  and  t5  =  sAt,  5  =  0, 1,  2, . . . ,  5.  The 

midpoint  between  two  mesh  points  in  x  can  be  interpolated  linearly  by  xr+i  =  The 

approximation  (4.21)  may  be  utilized  to  obtain 


Pr+1  ~  Pr 

At 


ps+\  _  J?s+l 
r+5  T~\ 

Ax 


(4.22) 


where  p^1  —  p*r  denotes  the  change  of  the  p.d.f.  p  w.r.t.  the  change  in  x  G  x  from  one  rnesh 
point  to  the  proceeding  mesh  point.  Letting  A  =  ct  if  x  G  C  and  A  =  7  if  x  G  e,  may 

r+  2 

be  defined  as 
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r+k 


f^s+l  T)-s+l  1 
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(4.23) 
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where  A) 1 1  = 


exp(Tr’  +  1)- 


T',+  1  = 


r+ 


A2 


To  preserve  the  zero-flux  boundary  conditions  of  F(x|xm*n  G  x,^+i)  —  0  and  F(x \xm(U.  G 

x,  ts+i)  =  07  let  =  0  and  Fs+l  =  0.  Equation  (4.22)  can  be  solved  using  the  system 
2  r_'"  2 


— (bs+lT)s+1  -4-  (hs+iT)s+l  —  fhs+1  T)s+1 
Vh,r  Pr+ 1  ^  Wo,r  Pr  V'-l.rPr- 1 


7^r 


(4.24) 
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where 


K1  = 


«+i  _ 


A £ 

Ax 
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4  Ax  r 
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Ax 
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4.2.10  Conclusions  to  Date 

The  results  found  in  [7]  were  successfully  duplicated  using  a  modified  version  of  the  model. 
The  coupled  SDE  model  proposed  in  [7]  can  be  converted  to  an  SDE  model  with  smooth 
coefficients.  Equation  (4.6)  orginally  possessed  smooth  coefficients,  and  to  change  (4.5) 
to  be  composed  of  only  smooth  coefficients  involved  scaling  the  difference'  in  two  agents’ 
characteristic  values  by  the  connectivity  between  those  two  agents  if  the  connectivity  is 
positive.  Other  methods  of  smoothing  the  coefficients  were  considered  but  did  not  fulfill  the 
requirement  of  homophily  between  agents.  Changing  the  system  (4.5)-(4.6)  to  (4.9)-(4.10) 
we  found  no  discernable  change  in  model  behavior  in  both  deterministic  and  stochastic  cases 
and  no  change  in  the  effect  of  parameters  ft  and  b  or  variables  a  and  7  on  the  system. 

With  smooth  coefficients,  the  SDE  fulfills  the  definition  of  a  diffusion  process,  thus  al¬ 
lowing  an  equivalent  multidimensional  Fokker-Planck  (FP)  model  to  be  formulated  from  the 
coupled  SDE  system  (4.9)  -  (4  10).  Usage  of  the  FP  model  is  nearly  ubiquitous  in  many 
scientific  fields,  and  current  research  is  being  performed  to  further  develop  and  understand 
the  FP  equation  [3].  The  FP  formulation  possesses  drift  coefficients  as  described  in  (4.20) 
and  diffusion  coefficients  equal  to  a  in  J\fm  dimensions  of  the  process  and  7  in  an  additional 
—  1)  dimensions.  It  is  possible  to  approximate  the  FP  model  by  a  Finite  Difference 
Scheme;  however,  the  ability  and  computational  efficiency  of  such  a  method,  while  still  under 
investigation,  are  less  than  promising.  In  response  to  the  conclusions  of  these  investigations, 
it  was  decided  that  additional  model  complexity  will  be  incorporated  and  delays  will  be 
added  directly  to  the  SDE  model. 


4.2.11  Social  Network  Model  with  Delay 

In  a  first  step  the  version  of  the  stochastic  social  network  model  discussed  above  and  in  [8] 
can  be  generalized  with  a  single  discrete  delay  is 


Ci(t)  =  — - -  ei{if ,t  -  T)[Ci’(t  -  t)  -Ci{t-r)]  (4.25) 

2-^i'£Ai(t-T)  eU  )  )  ,,6i4.(t_T) 
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where  /(£)  =  2e~b *  —  1  and  <fr  =  ,t  —  r).  This  system  is  for  i  =  where  m  is 

the  number  of  agents  considered,  C*(£)  denotes  the  value(s)  of  the  k  characteristics  of  the 
agents,  and  e(z,  z',f)  the  strength  of  the  connections  between  the  agents. 

We  are  currently  using  delay  equation  numerical  approximation  techniques  [1,  6]  devel¬ 
oped  under  previous  AFOSR  support  to  explore  features  of  these  models  with  delays.  We 
note  the  dimension  of  the  approximating  ODE  system  is  then  n  —  km  +  m2  =  m(k  +  to). 
Results  will  be  reported  in  [9]. 
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